function Threaded3DGammaTest()
%Threaded3DGammaTest
%
%This file will generate a few simple tests. 
%
%Created by:
% Jeff Craig
% January 2011

% Copyright 2011, Jeff Craig and Michael Jensen.
%  
% Threaded 3D Gamma is distributed under the terms of the Lesser GNU Public License. 
% 
%     This version of Threaded 3D Gamma is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
% 
% Threaded 3D Gamma is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
% 
% You should have received a copy of the GNU General Public License
% along with Threaded 3D Gamma.  If not, see
% <http://www.gnu.org/licenses/>.

%% Create variables
RefDose3D_org = ones(20,20,20);
EvalDose3D_org = ones(20,20,20);
RefDose2D_org = ones(20,20);
EvalDose2D_org = ones(20,20);
VoxelLength_org = 3; %mm
DoseDiff_org = 3; %percent
DTA_org = 3; %mm
nThreads_org = 4;
interp_org = 1; %no change

%% 3D - No change

[Gamma3D_same,gammax,gammay,gammaz] = Threaded3DGamma(RefDose3D_org, EvalDose3D_org, VoxelLength_org, DoseDiff_org, DTA_org, nThreads_org, interp_org);

%plot gamma
figure();imagesc(Gamma3D_same(:,:,10));

%% 2D - No change

[Gamma2D_same] = Threaded3DGamma(RefDose2D_org, EvalDose2D_org, VoxelLength_org, DoseDiff_org, DTA_org, nThreads_org, interp_org);

%plot gamma
figure();imagesc(Gamma2D_same(:,:));

%% 3D - Percent diff

%1 percent difference
RefDose = RefDose3D_org;
RefDose(5:15,5:15,5:15) = 1*100;

EvalDose = EvalDose3D_org;
EvalDose(5:15,5:15,5:15) = 1*99;
[Gamma_1percent,gammax,gammay,gammaz] = Threaded3DGamma(RefDose, EvalDose, VoxelLength_org, DoseDiff_org, DTA_org, nThreads_org, interp_org);

%plot gamma
figure();imagesc(Gamma_1percent(:,:,10));

%2.9 percent difference
RefDose = RefDose3D_org;
RefDose(5:15,5:15,5:15) = 1*100;

EvalDose = EvalDose3D_org;
EvalDose(5:15,5:15,5:15) = 1*98.1;
[Gamma_2_9percent] = Threaded3DGamma(RefDose, EvalDose, VoxelLength_org, DoseDiff_org, DTA_org, nThreads_org, interp_org);

%plot gamma
figure();imagesc(Gamma_2_9percent(:,:,10));

%3 percent difference
RefDose = RefDose3D_org;
RefDose(5:15,5:15,5:15) = 1*100;

EvalDose = EvalDose3D_org;
EvalDose(5:15,5:15,5:15) = 1*97;
[Gamma_3percent] = Threaded3DGamma(RefDose, EvalDose, VoxelLength_org, DoseDiff_org, DTA_org, nThreads_org, interp_org);

%plot gamma
figure();imagesc(Gamma_3percent(:,:,10));

%10 percent difference
RefDose = RefDose3D_org;
RefDose(5:15,5:15,5:15) = 1*100;

EvalDose = EvalDose3D_org;
EvalDose(5:15,5:15,5:15) = 1*10;
[Gamma_10percent] = Threaded3DGamma(RefDose, EvalDose, VoxelLength_org, DoseDiff_org, DTA_org, nThreads_org, interp_org);

%plot gamma
figure();imagesc(Gamma_10percent(:,:,10));

%% 2D - Percent diff

%1 percent difference
RefDose = RefDose2D_org;
RefDose(5:15,5:15) = 1*100;

EvalDose = EvalDose3D_org;
EvalDose(5:15,5:15) = 1*99;
[Gamma_1percent] = Threaded3DGamma(RefDose, EvalDose, VoxelLength_org, DoseDiff_org, DTA_org, nThreads_org, interp_org);

%plot gamma
figure();imagesc(Gamma_1percent(:,:));

%2.9 percent difference
RefDose = RefDose2D_org;
RefDose(5:15,5:15) = 1*100;

EvalDose = EvalDose2D_org;
EvalDose(5:15,5:15) = 1*98.1;
[Gamma_2_9percent] = Threaded3DGamma(RefDose, EvalDose, VoxelLength_org, DoseDiff_org, DTA_org, nThreads_org, interp_org);

%plot gamma
figure();imagesc(Gamma_2_9percent(:,:));

%3 percent difference
RefDose = RefDose2D_org;
RefDose(5:15,5:15) = 1*100;

EvalDose = EvalDose2D_org;
EvalDose(5:15,5:15) = 1*97;
[Gamma_3percent] = Threaded3DGamma(RefDose, EvalDose, VoxelLength_org, DoseDiff_org, DTA_org, nThreads_org, interp_org);

%plot gamma
figure();imagesc(Gamma_3percent(:,:));

%10 percent difference
RefDose = RefDose2D_org;
RefDose(5:15,5:15) = 1*100;

EvalDose = EvalDose3D_org;
EvalDose(5:15,5:15) = 1*10;
[Gamma_10percent] = Threaded3DGamma(RefDose, EvalDose, VoxelLength_org, DoseDiff_org, DTA_org, nThreads_org, interp_org);

%plot gamma
figure();imagesc(Gamma_10percent(:,:));

%% 3D - shift in 1 direction

%1 cm difference
RefDose = RefDose3D_org;
RefDose(5:15,5:15,5:15) = 1*100;

EvalDose = EvalDose3D_org;
EvalDose(6:16,5:15,5:15) = 1*100;

VoxelLength = 1;%cm

[Gamma_1cm,gammax,gammay,gammaz] = Threaded3DGamma(RefDose, EvalDose, VoxelLength, DoseDiff_org, DTA_org, nThreads_org, interp_org);

%plot gamma
figure();imagesc(Gamma_1cm(:,:,10));

%2 cm difference
RefDose = RefDose3D_org;
RefDose(5:15,5:15,5:15) = 1*100;

EvalDose = EvalDose3D_org;
EvalDose(7:17,5:15,5:15) = 1*100;

VoxelLength = 1;%cm

[Gamma_2cm,gammax,gammay,gammaz] = Threaded3DGamma(RefDose, EvalDose, VoxelLength, DoseDiff_org, DTA_org, nThreads_org, interp_org);

%plot gamma
figure();imagesc(Gamma_2cm(:,:,10));

%3 cm difference
RefDose = RefDose3D_org;
RefDose(5:15,5:15,5:15) = 1*100;

EvalDose = EvalDose3D_org;
EvalDose(8:18,5:15,5:15) = 1*100;

VoxelLength = 1;%cm

[Gamma_3cm] = Threaded3DGamma(RefDose, EvalDose, VoxelLength, DoseDiff_org, DTA_org, nThreads_org, interp_org);

%plot gamma
figure();imagesc(Gamma_3cm(:,:,10));

%4 cm difference
RefDose = RefDose3D_org;
RefDose(5:15,5:15,5:15) = 1*100;

EvalDose = EvalDose3D_org;
EvalDose(9:19,5:15,5:15) = 1*100;

VoxelLength = 1;%cm

[Gamma_4cm,gammax,gammay,gammaz] = Threaded3DGamma(RefDose, EvalDose, VoxelLength, DoseDiff_org, DTA_org, nThreads_org, interp_org);

%plot gamma
figure();imagesc(Gamma_4cm(:,:,10));

%% 2D - shift in 1 direction

%1 cm difference
RefDose = RefDose2D_org;
RefDose(5:15,5:15) = 1*100;

EvalDose = EvalDose2D_org;
EvalDose(6:16,5:15) = 1*100;

VoxelLength = 1;%cm

[Gamma_1cm] = Threaded3DGamma(RefDose, EvalDose, VoxelLength, DoseDiff_org, DTA_org, nThreads_org, interp_org);

%plot gamma
figure();imagesc(Gamma_1cm(:,:));

%2 cm difference
RefDose = RefDose2D_org;
RefDose(5:15,5:15) = 1*100;

EvalDose = EvalDose2D_org;
EvalDose(7:17,5:15) = 1*100;

VoxelLength = 1;%cm

[Gamma_2cm] = Threaded3DGamma(RefDose, EvalDose, VoxelLength, DoseDiff_org, DTA_org, nThreads_org, interp_org);

%plot gamma
figure();imagesc(Gamma_2cm(:,:));

%3 cm difference
RefDose = RefDose2D_org;
RefDose(5:15,5:15) = 1*100;

EvalDose = EvalDose2D_org;
EvalDose(8:18,5:15) = 1*100;

VoxelLength = 1;%cm

[Gamma_3cm] = Threaded3DGamma(RefDose, EvalDose, VoxelLength, DoseDiff_org, DTA_org, nThreads_org, interp_org);

%plot gamma
figure();imagesc(Gamma_3cm(:,:));

%4 cm difference
RefDose = RefDose2D_org;
RefDose(5:15,5:15) = 1*100;

EvalDose = EvalDose2D_org;
EvalDose(9:19,5:15) = 1*100;

VoxelLength = 1;%cm

[Gamma_4cm] = Threaded3DGamma(RefDose, EvalDose, VoxelLength, DoseDiff_org, DTA_org, nThreads_org, interp_org);

%plot gamma
figure();imagesc(Gamma_4cm(:,:));

%% 3D - shift in 2 directions

%1 cm difference
RefDose = RefDose_org;
RefDose(5:15,5:15,5:15) = 1*100;

EvalDose = EvalDose_org;
EvalDose(6:16,6:16,5:15) = 1*100;

VoxelLength = 1;%cm

[Gamma_1cm] = Threaded3DGamma(RefDose, EvalDose, VoxelLength, DoseDiff_org, DTA_org, nThreads_org, interp_org);

%plot gamma
figure();imagesc(Gamma_1cm(:,:,10));

%2 cm difference
RefDose = RefDose_org;
RefDose(5:15,5:15,5:15) = 1*100;

EvalDose = EvalDose_org;
EvalDose(7:17,7:17,5:15) = 1*100;

VoxelLength = 1;%cm

[Gamma_2cm] = Threaded3DGamma(RefDose, EvalDose, VoxelLength, DoseDiff_org, DTA_org, nThreads_org, interp_org);

%plot gamma
figure();imagesc(Gamma_2cm(:,:,10));

%3 cm difference
RefDose = RefDose_org;
RefDose(5:15,5:15,5:15) = 1*100;

EvalDose = EvalDose_org;
EvalDose(8:18,8:18,5:15) = 1*100;

VoxelLength = 1;%cm

[Gamma_3cm] = Threaded3DGamma(RefDose, EvalDose, VoxelLength, DoseDiff_org, DTA_org, nThreads_org, interp_org);

%plot gamma
figure();imagesc(Gamma_3cm(:,:,10));

%4 cm difference
RefDose = RefDose_org;
RefDose(5:15,5:15,5:15) = 1*100;

EvalDose = EvalDose_org;
EvalDose(9:19,9:19,5:15) = 1*100;

VoxelLength = 1;%cm

[Gamma_4cm] = Threaded3DGamma(RefDose, EvalDose, VoxelLength, DoseDiff_org, DTA_org, nThreads_org, interp_org);

%plot gamma
figure();imagesc(Gamma_4cm(:,:,10));

%% 2D - shift in 2 directions

%1 cm difference
RefDose = RefDose_org;
RefDose(5:15,5:15) = 1*100;

EvalDose = EvalDose_org;
EvalDose(6:16,6:16) = 1*100;

VoxelLength = 1;%cm

[Gamma_1cm] = Threaded3DGamma(RefDose, EvalDose, VoxelLength, DoseDiff_org, DTA_org, nThreads_org, interp_org);

%plot gamma
figure();imagesc(Gamma_1cm(:,:));

%2 cm difference
RefDose = RefDose_org;
RefDose(5:15,5:15) = 1*100;

EvalDose = EvalDose_org;
EvalDose(7:17,7:17) = 1*100;

VoxelLength = 1;%cm

[Gamma_2cm] = Threaded3DGamma(RefDose, EvalDose, VoxelLength, DoseDiff_org, DTA_org, nThreads_org, interp_org);

%plot gamma
figure();imagesc(Gamma_2cm(:,:));

%3 cm difference
RefDose = RefDose_org;
RefDose(5:15,5:15) = 1*100;

EvalDose = EvalDose_org;
EvalDose(8:18,8:18) = 1*100;

VoxelLength = 1;%cm

[Gamma_3cm] = Threaded3DGamma(RefDose, EvalDose, VoxelLength, DoseDiff_org, DTA_org, nThreads_org, interp_org);

%plot gamma
figure();imagesc(Gamma_3cm(:,:));

%4 cm difference
RefDose = RefDose_org;
RefDose(5:15,5:15) = 1*100;

EvalDose = EvalDose_org;
EvalDose(9:19,9:19) = 1*100;

VoxelLength = 1;%cm

[Gamma_4cm] = Threaded3DGamma(RefDose, EvalDose, VoxelLength, DoseDiff_org, DTA_org, nThreads_org, interp_org);

%plot gamma
figure();imagesc(Gamma_4cm(:,:));

%% 3D - shift in 3 directions

%1 cm difference
RefDose = RefDose_org;
RefDose(5:15,5:15,5:15) = 1*100;

EvalDose = EvalDose_org;
EvalDose(6:16,6:16,6:16) = 1*100;

VoxelLength = 1;%cm

[Gamma_1cm] = Threaded3DGamma(RefDose, EvalDose, VoxelLength, DoseDiff_org, DTA_org, nThreads_org, interp_org);

%plot gamma
figure();imagesc(Gamma_1cm(:,:,10));

%2 cm difference
RefDose = RefDose_org;
RefDose(5:15,5:15,5:15) = 1*100;

EvalDose = EvalDose_org;
EvalDose(7:17,7:17,7:17) = 1*100;

VoxelLength = 1;%cm

[Gamma_2cm] = Threaded3DGamma(RefDose, EvalDose, VoxelLength, DoseDiff_org, DTA_org, nThreads_org, interp_org);

%plot gamma
figure();imagesc(Gamma_2cm(:,:,10));

%3 cm difference
RefDose = RefDose_org;
RefDose(5:15,5:15,5:15) = 1*100;

EvalDose = EvalDose_org;
EvalDose(8:18,8:18,8:18) = 1*100;

VoxelLength = 1;%cm

[Gamma_3cm] = Threaded3DGamma(RefDose, EvalDose, VoxelLength, DoseDiff_org, DTA_org, nThreads_org, interp_org);

%plot gamma
figure();imagesc(Gamma_3cm(:,:,10));

%4 cm difference
RefDose = RefDose_org;
RefDose(5:15,5:15,5:15) = 1*100;

EvalDose = EvalDose_org;
EvalDose(9:19,9:19,9:19) = 1*100;

VoxelLength = 1;%cm

[Gamma_4cm] = Threaded3DGamma(RefDose, EvalDose, VoxelLength, DoseDiff_org, DTA_org, nThreads_org, interp_org);

%plot gamma
figure();imagesc(Gamma_4cm(:,:,10));
%% 2D shift and percent change

%1cm/1% difference
RefDose = RefDose_org;
RefDose(5:15,5:15,5:15) = 1*100;

EvalDose = EvalDose_org;
EvalDose(6:16,6:16,5:15) = 1*99;

VoxelLength = 1;%cm

[Gamma_1cm] = Threaded3DGamma(RefDose, EvalDose, VoxelLength, DoseDiff_org, DTA_org, nThreads_org, interp_org);

%plot gamma
figure();imagesc(Gamma_1cm(:,:,10));

%2cm/2% difference
RefDose = RefDose_org;
RefDose(5:15,5:15,5:15) = 1*100;

EvalDose = EvalDose_org;
EvalDose(7:17,7:17,5:15) = 1*98;

VoxelLength = 1;%cm

[Gamma_2cm] = Threaded3DGamma(RefDose, EvalDose, VoxelLength, DoseDiff_org, DTA_org, nThreads_org, interp_org);

%plot gamma
figure();imagesc(Gamma_2cm(:,:,10));

%3cm/3% difference
RefDose = RefDose_org;
RefDose(5:15,5:15,5:15) = 1*100;

EvalDose = EvalDose_org;
EvalDose(8:18,8:18,5:15) = 1*97;

VoxelLength = 1;%cm

[Gamma_3cm] = Threaded3DGamma(RefDose, EvalDose, VoxelLength, DoseDiff_org, DTA_org, nThreads_org, interp_org);

%plot gamma
figure();imagesc(Gamma_3cm(:,:,10));

%4cm/10% difference
RefDose = RefDose_org;
RefDose(5:15,5:15,5:15) = 1*90;

EvalDose = EvalDose_org;
EvalDose(9:19,9:19,5:15) = 1*100;

VoxelLength = 1;%cm

[Gamma_4cm] = Threaded3DGamma(RefDose, EvalDose, VoxelLength, DoseDiff_org, DTA_org, nThreads_org, interp_org);

%plot gamma
figure();imagesc(Gamma_4cm(:,:,10));


%% gamma with interp

%No interp
RefDose = RefDose_org;
RefDose(5:15,5:15,5:15) = 1*100;

EvalDose = EvalDose_org;
EvalDose(5:15,5:15,5:15) = 1*99;

interp = 1;
tic
[Gamma_100interp] = Threaded3DGamma(RefDose, EvalDose, VoxelLength_org, DoseDiff_org, DTA_org, nThreads_org, interp);
toc
%plot gamma
figure();imagesc(Gamma_100interp(:,:,10));


%interp to 50% of original size
RefDose = RefDose_org;
RefDose(5:15,5:15,5:15) = 1*100;

EvalDose = EvalDose_org;
EvalDose(5:15,5:15,5:15) = 1*99;

interp = .5;
tic
[Gamma_50interp] = Threaded3DGamma(RefDose, EvalDose, VoxelLength_org, DoseDiff_org, DTA_org, nThreads_org, interp);
toc
%plot gamma
figure();imagesc(Gamma_50interp(:,:,10));


%interp to 25% of original size
RefDose = RefDose_org;
RefDose(5:15,5:15,5:15) = 1*100;

EvalDose = EvalDose_org;
EvalDose(5:15,5:15,5:15) = 1*99;

interp = .25;
tic
[Gamma_25interp] = Threaded3DGamma(RefDose, EvalDose, VoxelLength_org, DoseDiff_org, DTA_org, nThreads_org, interp);
toc
%plot gamma
figure();imagesc(Gamma_25interp(:,:,10));


%interp to 5% of original size
RefDose = RefDose_org;
RefDose(5:15,5:15,5:15) = 1*100;

EvalDose = EvalDose_org;
EvalDose(5:15,5:15,5:15) = 1*99;

interp = .05;
tic
[Gamma_5interp] = Threaded3DGamma(RefDose, EvalDose, VoxelLength_org, DoseDiff_org, DTA_org, nThreads_org, interp);
toc
%plot gamma
figure();imagesc(Gamma_5interp(:,:,10));


%% timing gamma

avg = 0;

%1 thread
RefDose = RefDose_org;
RefDose(5:15,5:15,5:15) = 1*100;

EvalDose = EvalDose_org;
EvalDose(5:15,5:15,5:15) = 1*99;

nThreads = 1;
for i= 1:100
    tic
    [Gamma_1thread] = Threaded3DGamma(RefDose, EvalDose, VoxelLength_org, DoseDiff_org, DTA_org, nThreads, interp_org);
    avg = avg + toc;
end
avg_1 = avg / 100


%2 thread
RefDose = RefDose_org;
RefDose(5:15,5:15,5:15) = 1*100;

EvalDose = EvalDose_org;
EvalDose(5:15,5:15,5:15) = 1*99;

nThreads = 2;
for i= 1:100
    tic
    [Gamma_2thread] = Threaded3DGamma(RefDose, EvalDose, VoxelLength_org, DoseDiff_org, DTA_org, nThreads, interp_org);
    avg = avg + toc;
end
avg_2 = avg / 100


%4 thread
RefDose = RefDose_org;
RefDose(5:15,5:15,5:15) = 1*100;

EvalDose = EvalDose_org;
EvalDose(5:15,5:15,5:15) = 1*99;

nThreads = 4;
for i= 1:100
    tic
    [Gamma_4thread] = Threaded3DGamma(RefDose, EvalDose, VoxelLength_org, DoseDiff_org, DTA_org, nThreads, interp_org);
    avg = avg + toc;
end
avg_4 = avg / 100


%5 thread
RefDose = RefDose_org;
RefDose(5:15,5:15,5:15) = 1*100;

EvalDose = EvalDose_org;
EvalDose(5:15,5:15,5:15) = 1*99;

nThreads = 5;
for i= 1:100
    tic
    [Gamma_5thread] = Threaded3DGamma(RefDose, EvalDose, VoxelLength_org, DoseDiff_org, DTA_org, nThreads, interp_org);
    avg = avg + toc;
end
avg_5 = avg / 100


%8 thread
RefDose = RefDose_org;
RefDose(5:15,5:15,5:15) = 1*100;

EvalDose = EvalDose_org;
EvalDose(5:15,5:15,5:15) = 1*99;

nThreads = 8;
for i= 1:100
    tic
    [Gamma_8thread] = Threaded3DGamma(RefDose, EvalDose, VoxelLength_org, DoseDiff_org, DTA_org, nThreads, interp_org);
    avg = avg + toc;
end
avg_8 = avg / 100
